function [Nx,elem, node, int_elem, out_elem, immersed_inter, Dirichlet] = newest(h,Lx,By,Rx,Ty,r);
% function [Nx,elem, node, int_elem, out_elem, immersed_inter, Dirichlet] = newest(h,Lx,By,Rx,Ty,r)

if nargin < 6
    r  = 0.5;                %radius
end
if nargin <2
    Lx = -1;                 %left boundary
    By = -1;                 %bottom boundary
    Rx = 1 ;                  %right boundary     
    Ty = 1 ;                  %top boundary
end
if nargin < 1
    h = 1/20;
end
 
%%%%%%%%%%%%%%%%%%%%%%%
% Nx=199;
% Ny=199;
% hx=(Rx-Lx)/(Nx-1);
% hy=(Ty-Bx)/(Ny-1);
Nx = (Rx-Lx)/h+1;        %the number of node in  x direction
Ny = (Ty-By)/h+1;        %the number of node in  y direction
N  = Nx*Ny;              %the number of grid point
A  = [0,0];              % coordinate of the circle  center
%%%%%%%%%%the node coordinate in Cartesian mesh%%%%%%%%
   X  =  Lx:h:Rx;
   Y  =  By:h:Ty;
 node =  zeros(N,2);     %the number of grid points
for i = 1:Ny
    node((i-1)*Nx+1:i*Nx,1) = X(1:Nx);
    node((i-1)*Nx+1:i*Nx,2) = Y(i);
end
   %%%%%%%%%for check %%%
 flag_node = node;

%%%%%%%% refine grid%%%%%%%%%%%%%%%%%%%%
   %%%%%the interface point%%%%%%%
   %%%%%%%%row first then the column  %%%
x_min = ceil((A(1,1)-r-Lx)/h) ;      % left
x_max = floor((A(1,1)+r-Lx)/h);     % right

   %%%%%%%%correct the grid along axis x with y%%%%%%
for i = x_min:x_max
    y1 = sqrt(r^2-X(i+1)^2) ;
    y2 = -y1;
    m1 = ceil((y1-By)/(h/2));
    m2 = ceil((y2-By)/(h/2));
    
    if mod(m1,2) == 0
        
                node(m1/2*Nx+i+1,2) = y1;
           
       else node((m1-1)/2*Nx+i+1,2) = y1;
               
    end
   
    if mod(m2,2) == 0
                node(m2/2*Nx+i+1,2) =  y2;
    else    node((m2-1)/2*Nx+i+1,2) =  y2;
    end

%     
end
  
 y_min  = ceil((A(1,2)-r-By)/h);    % bottom
 y_max  = floor((A(1,2)+r-By)/h);   %top
 for i  = y_min:y_max
     x1 = sqrt(r^2-Y(i+1)^2);
     x2 = -x1;
     n1 = ceil(2*(x1-Lx)/h);
     n2 = ceil(2*(x2-Lx)/h);
     
     if mod(n1,2) == 0 
               node(i*Nx+n1/2+1,1) = x1;
        else node(i*Nx+(n1+1)/2,1) = x1;
     end
     
     if mod(n2,2) == 0 
               node(i*Nx+n2/2+1,1) = x2;
        else node(i*Nx+(n2+1)/2,1) = x2;
     end
 end
     
%%%%%%%%find the point  whos coordinate are changed in two
%%%%%%%%direction%%%%%%%%%

      flag = node - flag_node;
       C = flag(:,1) .* flag(:,2);
indices1 = find(C);
       I = size(indices1,1);
%%%%%%%%%%correct the grid point coordinate  use the nearer interface point%%%%%%%%
for  i = 1:I
    dx = abs(node(indices1(i),1)-flag_node(indices1(i),1));
    dy = abs(node(indices1(i),2)-flag_node(indices1(i),2));
    if dx>=dy
        node(indices1(i),1) =  flag_node(indices1(i),1)   ;
    else      
       node(indices1(i),2)  =  flag_node(indices1(i),2)   ;
    end
end


%%%%%%triangle elememt %%%%%%%
Ntri=2*(Nx-1)*(Ny-1);                %the number of triangle
elem=zeros(Ntri,3);                  %node script of every triangle
for i=1:Ny-1
             elem((i-1)*(Nx-1)+1:i*(Nx-1),1)  =    (i-1)*Nx+2:i*Nx ;
             elem((i-1)*(Nx-1)+1:i*(Nx-1),2)  =    i*Nx+2:(i+1)*Nx ;
             elem((i-1)*(Nx-1)+1:i*(Nx-1),3)  =  (i-1)*Nx+1:i*Nx-1 ;
    elem((i+Nx-2)*(Nx-1)+1:(i+Nx-1)*(Nx-1),1) =  i*Nx+1:(i+1)*Nx-1 ; 
    elem((i+Nx-2)*(Nx-1)+1:(i+Nx-1)*(Nx-1),2) =  (i-1)*Nx+1:i*Nx-1 ;
    elem((i+Nx-2)*(Nx-1)+1:(i+Nx-1)*(Nx-1),3) =    i*Nx+2:(i+1)*Nx ; 
end

%%%%%%%find the point  whos coordinate are changed %%%%%%%%%



%%%%%%%to decide which point is on interface%%%%%%
   nup   =   (Nx-1)*(Ny-1)                  ;
   absx  =   abs(node(:,1))                 ;
   absy  =   abs(node(:,2))                 ;
   dir_x =   abs(absx-sqrt(r^2-absy.*absy)) ;
   dir_y =   abs(absy-sqrt(r^2-absx.*absx)) ;
 epsilon =   1.0E-10                        ;
   xandy = (dir_x<epsilon)+(dir_y<epsilon)  ;
indices2 = find(xandy)                      ;

% equalx=find(absx==sqrt(r^2-absy.*absy));
% equaly=find(absy==sqrt(r^2-absx.*absx));
% equalsy=find(absy==sqrt(r^2-absx.*absx)&&absy==sqrt(r^2-absx.*absx));
% 
% d=r^2*ones(size(a,1),1);
% indices2=find(d==c);
   nc = size(indices2,1);
for i = 1:nc
        j = indices2(i) ;
      row = floor(j/Nx) ;
    colmn =   mod(j,Nx) ;
    %%%%%%%%left  bottom %%%%%%%%
     if xandy(j-Nx-1,1) == 0||(xandy(j-1,1)+xandy(j-Nx,1)>=1)%(flag_node(j-Nx-1)==node(j-Nx-1))*(flag_node(j-Nx-1,2)==node(j-Nx-1,2))
        if  (node(j   ,1) - node(j-Nx-1,1))^2 + (node(j   ,2) - node(j-Nx-1,2))^2>...
            (node(j-Nx,1) - node(j-1   ,1))^2 + (node(j-Nx,2) - node(j-1   ,2))^2 %||xandy(j-1,1)&&xandy(j-Nx,1)%(node(j-1,1)^2+node(j-1,2)^2==r^2)*(node(j-Nx,1)^2+node(j-Nx,2)^2==r^2)
            elem(j-Nx-row    ,1) = j-Nx-1;
            elem(j-Nx-row    ,2) = j-Nx  ;
            elem(j-Nx-row    ,3) = j-1   ;
            elem(j-Nx-row+nup,1) = j     ;
            elem(j-Nx-row+nup,2) = j-1  ;
            elem(j-Nx-row+nup,3) = j-Nx   ;
        end
    end
    
 %%%%%%%%%%%%%%right  bottom %%%%%%%%
      
    if  (node(j     ,1) - node(j-Nx+1,1))^2 + (node(j   ,2) - node(j-Nx+1,2))^2<...
        (node(j-Nx  ,1) - node(j+1   ,1))^2 + (node(j-Nx,1) - node(j+1   ,1))^2||...
        (xandy(j-Nx+1,1)&&xandy(j+1,1) == 0&&xandy(j-Nx,1) == 0)%flag_node(j-Nx+1)~=node(j-Nx+1)%(node(j,1)^2+node(j,2)^2==r^2)*(node(j-Nx+1,1)^2+node(j-Nx+1,2)^2==r^2)
         elem(j-Nx-row+1     ,1) = j-Nx   ;
         elem(j-Nx-row+1     ,2) = j-Nx+1 ;
         elem(j-Nx-row+1     ,3) = j      ;
         elem(j-Nx-row+nup+1 ,1) = j+1    ;
         elem(j-Nx-row+nup+1 ,2) = j      ;
         elem(j-Nx-row+nup+1 ,3) = j+1-Nx ;
    end
%    end
    
    
 
     
      %%%%%%%%left  top %%%%%%%%
    if  (node(j     ,1) - node(j+Nx-1,1))^2 + (node(j   ,2) - node(j+Nx-1,2))^2<...
        (node(j+Nx  ,1) - node(j-1   ,1))^2 + (node(j+Nx,2) - node(j-1   ,2))^2||...
         (xandy(j+Nx+1,1)&&xandy(j-1,1) == 0&&xandy(j+Nx,1) == 0)%(node(j,1)^2+node(j,2)^2==r^2)*(node(j+Nx-1,1)^2+node(j+Nx-1,2)^2==r^2)
        elem(j-row-1     ,1) = j-1    ;
        elem(j-row-1     ,2) = j      ;
        elem(j-row-1     ,3) = j-1+Nx ;
        elem(j-row-1+nup ,1) = j+Nx   ;
        elem(j-row-1+nup ,2) = j+Nx-1 ;
        elem(j-row-1+nup ,3) = j      ;
    end
   % end
 
%     
%     
%     
    
      %%%%%%%%right  top %%%%%%%%
    if xandy(j+Nx+1,1) == 0||(xandy(j+1,1)+xandy(j+Nx,1)>=1)%(flag_node(j+Nx+1)==node(j+Nx+1))*(flag_node(j+Nx+1,2)==node(j+Nx+1,2))
       if (node(j   ,1) - node(j+Nx+1,1))^2 + (node(j   ,2)  -node(j+Nx+1,2))^2>...
          (node(j+Nx,1) - node(j+1   ,1))^2 + (node(j+Nx,2) - node(j+1   ,2))^2%||(node(j+1,1)^2+node(j+1,2)^2==r^2)*(node(j+Nx,1)^2+node(j+Nx,2)^2==r^2)
           elem(j-row    ,1) = j      ;
           elem(j-row    ,2) = j+1    ;
           elem(j-row    ,3) = j+Nx   ;
           elem(j-row+nup,1) = j+Nx+1 ;
           elem(j-row+nup,2) = j+Nx   ;
           elem(j-row+nup,3) = j+1    ;
       end
    end

end


%%%%%%%%%%%%%%interface%%%%%%%%%%%%%%%%%%%%%
     immersed_inter = zeros(nc,2)       ;
immersed_inter(1,1) = mod(indices2(1),N);

 for i = 2:nc
     j = immersed_inter(i-1,1);

    if node(j,1) >= 0 && node(j,1) <= r
        if node(j,2) >= (-r) && node(j,2) < 0
            if  xandy(j+1,1)%node(j+1,1)^2+node(j+1,2)^2==r^2  %%%%%%%% right %%%%%%
                   immersed_inter(i,1) =  j+  1 ;
            
            elseif  xandy(j+Nx+1)%node(j+1+Nx,1)^2+node(j+1+Nx,2)^2==r^2 %%%%%%%% right top %%%%%%
                   immersed_inter(i,1) = j+Nx+1 ;
            else   immersed_inter(i,1) = j+Nx   ;         %%%%%%%% top %%%%%%
             end
          
            
        elseif  node(j,2) >=0 && node(j,2) <= r
            if  xandy(j+Nx,1)%node(j+Nx,1)^2+node(j+Nx,2)^2==r^2  %%top
                   immersed_inter(i,1) = j+Nx   ;
            elseif   xandy(j+Nx-1,1)%node(j-1+Nx,1)^2+node(j-1+Nx,2)^2==r^2 %%%%%%%% left top %%%%%%
                   immersed_inter(i,1) = j+Nx-1 ;
            else   immersed_inter(i,1) = j-1    ;             %%%%%left

            end 
        
        end 
    end 
            
        
   if node(j,1) >= (-r) && node(j,1) < 0
       if node(j,2) >=0 && node(j,2) <= r 
           if  xandy(j-1,1)%node(j-1,1)^2+node(j-1,2)^2==r^2  %%%%%%%% left %%%%%%
                     immersed_inter(i,1) = j-1    ;
             
           elseif  xandy(j-Nx-1,1)%node(j-1-Nx,1)^2+node(j-1-Nx,2)^2==r^2 %%%%%%%% left bottom %%%%%%
                     immersed_inter(i,1) = j-Nx-1 ;
                    
           else      immersed_inter(i,1) = j-Nx   ;       %%%%%%%% bottom %%%%%%
               
           end
            
           elseif node(j,2)>=(-r)&&node(j,2)<0
                if xandy(j-Nx,1)%node(j-Nx,1)^2+node(j-Nx,2)^2==r^2  %%bottom
                     immersed_inter(i,1) = j-Nx   ;
              
             elseif  xandy(j-Nx+1,1)%node(j+1-Nx,1)^2+node(j+1-Nx,2)^2==r^2 %%%%%%%% right bottom %%%%%%
                     immersed_inter(i,1) = j-Nx+1 ; 
                else immersed_inter(i,1) = j+1    ;            %%%%%right

                end 
        
       end 
   end 

     immersed_inter(i-1,2) = immersed_inter(i,1);
 end 
immersed_inter(nc,2)=immersed_inter(1,1);



%%%%%%%%%the out boundary%%%%%%%%%
Dirichlet=zeros(2*(Nx+Ny-2),2);
%%%%%%%left%%%%%%%%
Dirichlet(1:Ny-1,1)=1:Nx:(Ny-2)*Nx+1;
Dirichlet(1:Ny-1,2)=Nx+1:Nx:(Ny-1)*Nx+1;
%%%%%%%%%%%top%%%%%%%%%%
Dirichlet(Ny:Ny+Nx-2,1)=(Ny-1)*Nx+1:Nx*Ny-1;
Dirichlet(Ny:Ny+Nx-2,2)=(Ny-1)*Nx+2:Nx*Ny;
%%%%%%%%%%right%%%%%%%%%%%
Dirichlet(Ny+Nx-1:2*Ny+Nx-3,1)=Nx*Ny:-Nx:2*Nx;
Dirichlet(Ny+Nx-1:2*Ny+Nx-3,2)=Nx*(Ny-1):-Nx:Nx;
%%%%%%%%%%%bottom%%%%%%%%%%
Dirichlet(2*Ny+Nx-2:2*(Nx+Ny-2),1)=Nx:-1:2;
Dirichlet(2*Ny+Nx-2:2*(Nx+Ny-2),2)=Nx-1:-1:1;


 
j = 1;
k = 1;
out_elem = zeros(1,3);
int_elem = zeros(1,3);
for i=1:Ntri
    if node(elem(i,1),1)^2+node(elem(i,1),2)^2+...
       +node(elem(i,2),1)^2+node(elem(i,2),2)^2+...
       +node(elem(i,3),1)^2+node(elem(i,3),2)^2>3*r^2
         out_elem(j,:) = elem(i,:) ;
                     j = j+1       ;
    else int_elem(k,:) = elem(i,:) ;
                     k = k+1       ;
    end
end
% elem = label(node,elem);
% figure(1)
% trisurf(elem,node(:,1),node(:,2),ones(size(node,1),1))
% view(2), axis equal, axis on,  title('triangle mesh') , drawnow;

% figure(2)
% plot(node(immersed_inter(:,1),1),node(immersed_inter(:,1),2),'-*');
%  title('the immersed interface point ')
% figure(3)
% %out_elem = label(node,out_elem);
% trisurf(out_elem,node(:,1),node(:,2),ones(size(node,1),1))
% view(2), axis equal, axis on, title('outer triangle element mesh') ,  drawnow;
%  
% figure(4)
% trisurf(int_elem,node(:,1),node(:,2),ones(size(node,1),1))
% view(2), axis equal, axis on,  title('inner triangle element mesh'),  drawnow;